library(dplyr)
library(knitr)
library(maptools)
library(rgdal)
library(TreeSegmentation)
library(sp)
library(ggplot2)
library(rgl)
library(clue)
library(lidR)
knit_hooks$set(webgl = hook_webgl)
opts_chunk$set(warning=F,message=F)

#set color ramp for treeID
col = pastel.colors(200)

#set data paths
path_to_tiles="../data/2017/Lidar/cropped_"

#set cores
cores<-4
#cores<-15

1 Load in ground-truth

shps<-list.files("../data/ITCs/test/",pattern=".shp",full.names = T)
itcs<-lapply(shps,readOGR,verbose=F)

names(itcs)<-sapply(itcs,function(x){
  id<-unique(x$Plot_ID)
  return(id)
  })

2 Example Pipeline

2.1 Read in Data

ground_truth<-itcs[[2]]
fname<-get_tile_filname(ground_truth)
tile<-readLAS(paste("../data/2017/Lidar/cropped_",fname,sep=""))
tile@crs<-CRS("+init=epsg:32617")
plot(tile)

You must enable Javascript to view this page properly.

2.2 Confirm overlap

plot(extent(tile),col='red')
plot(extent(ground_truth),col='blue',add=T)

2.3 Compute Segmentation

silva<-silva2016(tile=tile,output="all")
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##  28.375   0.529  29.012 
## [1] "Creating tree polygons"
dalponte<-dalponte2016(tile=tile,output="all")
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##  31.769   0.407  32.308 
## [1] "Creating tree polygons"
li<-li2012(tile=tile,output="all")
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##   0.854   0.016   0.873 
## [1] "Creating tree polygons"
watershed_result<-watershed(tile=tile,output="all")
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
##    user  system elapsed 
##  39.756   0.449  40.367 
## [1] "Creating tree polygons"

2.4 View 3d segmentation

plot(silva$tile,color="treeID",col=col)

You must enable Javascript to view this page properly.

2.4.1 Overlay ground truth and predictions

plot(silva$convex)

plot(ground_truth,col='red')
plot(silva$convex,add=T)

#plot(dalponte$convex,add=T)

2.4.2 CHM versus predicted polygons

chm=canopy_model(silva$tile)
plot(chm,ext=extent(ground_truth))
plot(ground_truth,add=T,col='red')
plot(silva$convex,add=T)

Okay that’s not great, but let’s keep going for the moment.

2.4.3 Compare Methods

Silva v Dalponte

plot(silva$convex)
plot(dalponte$convex,add=T,col=rgb(0,0,255,20,maxColorValue=255))

plot(silva$convex)
plot(dalponte$convex,add=T,col=rgb(0,0,255,20,maxColorValue=255))
plot(li$convex,add=T,col=rgb(0,0,255,20,maxColorValue=255))

Li versus watershed

plot(li$convex)
plot(watershed_result$convex,add=T,col=rgb(0,0,255,20,maxColorValue=255))

2.5 Compute consensus

ptlist<-list(silva=silva$tile,li=li$tile,dalponte=dalponte$tile,watershed=watershed_result$tile)
system.time(consensus_result<-consensus(ptlist=ptlist,method="majority"))
##    user  system elapsed 
## 262.229  46.596 310.249
consensus_polygons<-get_convex_hulls(consensus_result,consensus_result@data$treeID)
plot(consensus_polygons,col=rgb(0,0,255,20,maxColorValue=255))

paste(length(consensus_polygons),"consensus clusters found")
## [1] "1244 consensus clusters found"

2.5.1 View consensus segmentation

#plot(consensus_result,color="treeID",col=col)

You must enable Javascript to view this page properly.

2.5.2 Compare Consensus

plot(silva$convex)
plot(dalponte$convex)
plot(li$convex)
plot(consensus_polygons,add=T,col=rgb(0,0,255,40,maxColorValue=255))

How many tree predictions?

library(pander)
unique_total<-sapply(c(ptlist,consensus_result),function(x) length(unique(x@data$treeID)))
df<-data.frame(Algorthm=c(names(ptlist),"consensus"),Total_Trees=as.numeric(unique_total))
pandoc.table(df,style="rmarkdown")
## 
## 
## | Algorthm  | Total_Trees |
## |:---------:|:-----------:|
## |   silva   |    1265     |
## |    li     |    1814     |
## | dalponte  |    1265     |
## | watershed |     944     |
## | consensus |    1245     |

2.6 Assign Trees

Each tree is assigned based on the maximum overlap. Pairwise membership is done using a Hungarian Algorithm. See clue::solve_LSAP.

assignment<-assign_trees(ground_truth,prediction=silva$convex)

2.7 Calculate Intersection over union

#loop through assignments and get jaccard statistic for each assignment pair
statdf<-calc_jaccard(assignment=assignment,ground_truth = ground_truth,prediction=silva$convex)
ggplot(statdf) + geom_histogram(aes(IoU)) + labs(x="Intersection over union") + theme_bw()

mean(statdf$IoU)
## [1] 0.2652395
median(statdf$IoU)
## [1] 0.2592368

3 As a wrapper for one tile, multiple methods, and consensus

results<-evaluate(ground_truth=itcs[[1]],algorithm = c("silva","dalponte","li"),path_to_tiles=path_to_tiles,compute_consensus = F,plot_results=T)

if(!is.null(results)){
  ggplot(results,aes(y=IoU,x=Method)) + geom_boxplot() + theme_bw()
results %>% group_by(Method) %>% summarize(mean=mean(IoU),median=median(IoU))
}

4 Parallel wrapper across all tiles

system.time(results_all<-evaluate_all(itcs=itcs[1:2],algorithm = c("dalponte","li","watershed"),path_to_tiles=path_to_tiles,cores=cores,extra=F,compute_consensus=F))
## [1] "ITC  has no overlap with cropped tile"
##    user  system elapsed 
##   0.121   0.009 216.312
#plot results
ggplot(results_all,aes(y=IoU,x=Method)) + geom_boxplot() + theme_bw()

#Compute test statistics
test_statistic<-results_all %>% group_by(Method) %>% summarize(mean=mean(IoU),median=median(IoU)) %>% mutate(Date=format(Sys.time(), "%m/%d/%y %H:%M:%S"))
test_statistic
## # A tibble: 3 x 4
##   Method     mean median Date             
##   <chr>     <dbl>  <dbl> <chr>            
## 1 dalponte  0.295  0.347 05/01/18 17:07:47
## 2 li        0.313  0.311 05/01/18 17:07:47
## 3 watershed 0.233  0.226 05/01/18 17:07:47
write.csv(test_statistic,"Results/results.csv",append = T,col.names = NA)